
clear
do "...\First.do"

********************************************************************************
* This do file uses patients who died prior to clinic closure to test for potentiel
* pre trend in the mortality outcome
* The test is done in several steps. 

** First, I identify all patient in the relevant clinics dying within 5 years of clinic closure
** Second, I predict new PCP (post closure if they were alive) from the surviving patients in the clinics
** Third, test if the predict new PCP SES explains mortality patterns before clinic closure
** This file does it on the individual level

** The first step is done in 8_pre_mortality_cliniclevel.do
********************************************************************************



********************************************************************************
* Load dataset
********************************************************************************

clear 
use "$work\temp_all.dta"

keep pnr year timing2 death ydernr low_ses mean_age mean_male mean_dk solo ku au sdu other N_doctors male age non_dk married
keep if inrange(timing2,-5,0)

g pre_sample=1


********************************************************************************
* Second, predict new PCP (post closure) for dying patients from the surviving patients in the clinics using the other patient. The prediction is done on the individal level
********************************************************************************

********************************************************************************
** Append with peer patients who are alive - they are used for the prediction
append using "$work\analysis_sample.dta"
********************************************************************************


********************************************************************************
* Timing variable
********************************************************************************
replace timing2=timing2 // Timing for pre-mortality patients
replace timing2=timing if timing!=. // Timing for non-pre-mortality patients
tab timing2 death, row

cap drop t2
g t2=timing2+4
tab timing2 t2
tab timing2 death if ydernr!=""



********************************************************************************
* GP identifier
********************************************************************************
g gp_fe2=gp_fe
replace gp_fe2=ydernr if pre_sample==1
drop ydernr 
rename gp_fe2 ydernr



********************************************************************************
* Keep relevant population
********************************************************************************
keep if inrange(age,40,70)
drop if year>2019
drop if min_year==2017 & pre_sample==0


********************************************************************************
** Keep a balanced pre-period for surviving patients
********************************************************************************
cap drop help N_pre_2
g help=1 if inrange(timing2,-4,0)
bys pnr: egen N_pre_2=total(help)
drop if N_pre_2!=5 & pre_sample!=1 



********************************************************************************
** Background characteristics - these are not yet defined for the dying population
********************************************************************************
replace non_dk=1 if (dk==0) & non_dk==.
replace non_dk=0 if (dk==1) & non_dk==.
tab non_dk

g _age0=age if timing2==-4
bys pnr: egen age0=max(_age0)
sum age0

********************************************************************************
* Number of patients
********************************************************************************

bys gp_fe2 timing2 low_ses: g N_patients=_N


********************************************************************************
* Predicting the chance of getting a low-SES PCP post closure for dying population
*********************************************************************************

cap drop res_low
areg help3 male i.age married low_ses non_dk if timing2==-4 & low_ses==1 & age<64 & N_patients>125, a(gp_fe2)
predict res_low if timing2==-4 & low_ses==1  & N_patients>125
sum res_low, d


cap drop res_high
areg help3 male i.age married low_ses non_dk if timing2==-4 & low_ses==0 & age<64 & N_patients>125, a(gp_fe2)
predict res_high if timing2==-4 & low_ses==0  & N_patients>125
sum res_high, d

corr help3 res_low if timing2==-4 & low_ses==1 & help3!=.
corr help3 res_high if timing2==-4 & low_ses==0 & help3!=.

*********************************************************************************
* Creates predicted SES quartiles
*********************************************************************************

cap drop _ses_low ses_low
xtile _ses_low=res_low, n(4)
bys pnr: egen ses_low=max(_ses_low)
tab ses_low
tab ses_low death

cap drop _ses_high ses_high
xtile _ses_high=res_high, n(4)
bys pnr: egen ses_high=max(_ses_high)
tab ses_high
tab ses_high death

tab ses_low if low_ses==1, m
tab ses_high if low_ses==0, m


forvalues v=1/4 {


cap drop SES_`v'
g SES_`v'=(ses_low>=`v') if ses_low!=. & low_ses==1
replace SES_`v'=0 if (ses_low<`v') & ses_low!=. & low_ses==1

replace SES_`v'=1 if (ses_high>=`v') & ses_high!=. & low_ses==0
replace SES_`v'=0 if (ses_high<`v') & ses_high!=. & low_ses==0


}


********************************************************************************
* 3. test if the predict new PCP SES explains mortality patterns before clinic closure
** Focus on the pre-period [-4;0]
** Cap the age at time -4 to 64, as this matches the patients in the surviving sample used in the analysis
********************************************************************************

cap postclose john
postfile john T timing ses beta upper lower using "$table\PreDeath_2", replace


forvalues v=2/4 {
areg death ib4.t2##SES_`v' male i.age married non_dk if low_ses==0 & inrange(t2,0,4) & age0<64, cluster(gp_fe2) a(gp_fe2)
mat A=r(table)
mat list A

post john (`v') (0) (0) (A[1,9]) (A[6,9]) (A[5,9])
post john (`v') (1) (0) (A[1,11]) (A[6,11]) (A[5,11])
post john (`v') (2) (0) (A[1,13]) (A[6,13]) (A[5,13])
post john (`v') (3) (0) (A[1,15]) (A[6,15]) (A[5,15])
post john (`v') (4) (0) (0) (.) (.)


areg death ib4.t2##SES_`v' male i.age married non_dk if low_ses==1 & inrange(t2,0,4) & age0<64, cluster(gp_fe2) a(gp_fe2)
mat A=r(table)
mat list A

post john (`v') (0) (1) (A[1,9]) (A[6,9]) (A[5,9])
post john (`v') (1) (1) (A[1,11]) (A[6,11]) (A[5,11])
post john (`v') (2) (1) (A[1,13]) (A[6,13]) (A[5,13])
post john (`v') (3) (1) (A[1,15]) (A[6,15]) (A[5,15])
post john (`v') (4) (1) (0) (.) (.)


areg death ib4.t2##SES_`v'##low_ses male i.age married non_dk if inrange(t2,0,4) & age0<64, cluster(gp_fe2) a(gp_fe2)
mat A=r(table)
mat list A

post john (`v') (0) (3) (A[1,37]) (A[6,37]) (A[5,37])
post john (`v') (1) (3) (A[1,41]) (A[6,41]) (A[5,41])
post john (`v') (2) (3) (A[1,45]) (A[6,45]) (A[5,45])
post john (`v') (3) (3) (A[1,49]) (A[6,49]) (A[5,49])
post john (`v') (4) (3) (0) (.) (.)

}

postclose john



preserve

clear 
use "$table\PreDeath_2"


replace timing=timing-4
replace beta=beta*100
replace upper=upper*100
replace lower=lower*100

replace upper=0 if timing==0
replace lower=0 if timing==0



tw 	(line beta timing if ses==1 & T==3, color(sienna%60)) ///
	(line beta timing if ses==1 & T==4, color(sienna)) ///
	(line beta timing if ses==0 & T==3, color(navy%60)) ///
	(line beta timing if ses==0 & T==4, color(navy) ///
	legend(order(1 "Low SES Q>=3" 2 "Low SES Q>=4" ///
		3 "High SES Q>=3" 4 "High SES Q>=4")  region(style(none)) rows(2) pos(6)) ///
	yscale(range(-0.4(0.2)0.4)) ylabel(-0.4(0.2)0.4) ///
	graphregion(color(white) ilcolor(white) lcolor(white)) bgcolor(white) ///
	ytitle("Percentage points") xtitle("Years since clinic closure")	)
graph export "$fig\Death_placebo_v2_idlevel.png", replace 
graph export "$fig\Death_placebo_v2_idlevel.pdf", replace 


tw	(rcap lower upper timing if ses==1 & T==4, color(sienna%50) lc(*0.5)) ///
	(line beta timing if ses==1 & T==4, color(sienna)) ///
	(rcap lower upper timing if ses==0 & T==4, color(navy%50) lc(*0.5)) ///
	(line beta timing if ses==0 & T==4, color(navy) ///
	legend(order(2 "Low SES" 4 "High SES") region(style(none)) rows(1) pos(6)) ///
	yscale(range(-0.4(0.2)0.4)) ylabel(-0.4(0.2)0.4) ///
	graphregion(color(white) ilcolor(white) lcolor(white)) bgcolor(white))
graph export "$fig\Death_placebo_Q5_v2_idlevel.png", replace 
graph export "$fig\Death_placebo_Q5_v2_idlevel.pdf", replace 

restore





